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Abstract 


The primary objective of this thesis is to investigate methods to improve the ability 
to maintain the inventory of orbital elements of Earth satellites during periods of extreme 
atmospheric disturbance brought on by severe solar activity. Existing tracking techniques 
do not account for such atmospheric dynamics. This can result in tracking errors of 
several seconds in predicted crossing time during periods of high geomagnetic activity for 
certain satellites. The reduction of these tracking errors is the principal goal of this thesis. 
Two techniques are examined for this purpose. In the first approach, density predicted 
from various atmospheric models is fit to the orbital decay rate for a number of satellites 
using a least-squares method. An orbital decay model is then developed that potentially 
could be used to reduce tracking errors by accounting for atmospheric changes. The 
second approach utilizes a Kalman filter to estimate the orbital decay rate of a satellite 
after every observation. The new information is then used to predict the next observation. 

Results from the first approach demonstrated the feasibility of building an orbital 
decay model based on predicted atmospheric density, which then could potentially be used 
to reduce tracking errors. Correlation of atmospheric density to orbital decay was as high 
as 0.88. However, it is clear that contemporary atmospheric models need further 
improvement in modeling density perturbations brought on by solar activity in the polar 
region. The second approach of Kalman filtering satellite orbital decay resulted in a 
dramatic reduction in tracking errors for certain satellites during severe solar storms. For 
example, in the limited cases studied, the reduction in tracking errors ranged from 79 to 25 
percent. 
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1. Introduction 


The tracking of artificial earth satellites began with the launch of Sputnik in 1957. 
During the Cold War it was imperative to distinguish existing satellite orbits from new 
launches to determine whether or not those new satellites posed a threat. Today, the 
environment is considerably different but the need still exists to be able to identify 
satellites. In addition, the U.S. has treaty agreements to predict satellite re-entry and 
notify the appropriate country where landfall is expected to occur. 1 Currently there are 
over 8,000 satellites in orbit, which is considerably more than even a few years ago. 
Consequently, congestion and collision avoidance are growing concerns and hence one of 
the continuing needs for space surveillance. An inventory of orbital information on every 
satellite is continuously maintained and hereinafter referred to as the Earth satellite 
catalog. 

Satellite tracking is achieved through the use of various radar installations such as 
the fence that is operated by the Naval Space Command (NSC) from Dahlgren, VA. 2 As a 
satellite passes through the fence, receivers detect a reflected signal and a crossing time is 
associated with the peak amplitude of the reflected signal. A predicted crossing time for 
the satellite being tracked is calculated using orbital theory. This predicted time is 
compared to the actual crossing time and if the difference is within a specified tolerance 
the satellite is considered identified. Otherwise, the satellite becomes an uncorrelated 
target (UCT). UCT’s can occur for a number of reasons including the launching of new 
objects, collisions between existing ones, explosions, propulsive maneuvers, 
approximations made in the geopotential of the Earth, or from disturbances in the 
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atmosphere resulting in a change in density and therefore a change in atmospheric drag. 
Increased drag decreases the orbital energy of a satellite and consequently the semi-major 
axis of the orbit. Since the square of the period of a satellite is proportional to the cube of 
the semi-major axis, a decrease in the semi-major axis will result in a decrease in the 

period of a satellite. Therefore, an increase in drag will cause a satellite to arrive earlier at 
the fence than predicted. 

Drag effects are incorporated in the orbit model by carrying, in addition to the six 

orbital elements, a seventh term, the rate of change in mean motion, or n. This term, also 
known as the orbital decay parameter (ODP), is used to absorb non-central forces such as 
drag, solar radiation pressure, forces due to higher order harmonics, thrusting, etc., that 
are not incorporated into the orbital model. 

A catalog entry for a new satellite is created shortly after launch when elements are 
initially calculated using orbit determination (OD) methods. Periodic updates to the 
orbital elements are required to maintain the catalog. This is accomplished by using a 
differential correction (DC) process of the orbital elements. This process requires many 
observations, spanning days, whereby residuals are created by differencing the actual 
observation from those predicted by the orbit model. In addition, DC of elements requires 
the solution of a set of j simultaneous linear equations in seven unknowns, where j is the 
number of observations. Since the number of observations is typically larger than the 
number of unknowns, there are more equations than unknowns and therefore no unique 
solution. To find the best solution, a batch least-squares (LS) process is employed. 
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providing the corrections to the elements that minimizes the square root of the arithmetic 
mean of the squared residuals, or root mean square (RMS). 3 

For low Earth orbiting (LEO) satellites, ODP is primarily an indication of 
atmospheric drag a satellite experiences. Because ODP is part of the LS solution to the 
DC process that can involve several days worth of observations, ODP will be an average 
indication of the drag acting on a satellite during the period of observation. Because 
atmospheric disturbances can have time spans on the order of hours, using the smoothed 
ODP to make fence crossing predictions has led to fence crossing errors greater than 10 
seconds. Such a large error is well outside the crossing tolerance of 2 seconds and 
therefore not conducive to maintaining the satellite catalog. 

The present research centers around reducing UCT’s that are caused by the 
atmospheric perturbations mentioned above. This is accomplished by studying several 
atmospheric models and, in the first approach, developing a simple drag or ODP model 
that could be applied to the current orbital equations and potentially reduce the number of 
UCT’s occurring. As a second approach, a Kalman filter is used to continuously update 
ODP after each observation. This updated ODP is then used in place of the batch ODP to 
reduce UCT’s. 

This thesis begins with a review of the governing equations, followed by 
descriptions of the atmospheric models studied. Next, a description of the radar fence and 
the data it generates is given. Thereafter the ODP model is presented, followed by the 
Kalman filter approach and concluding remarks. 
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Additional background information is provided in the appendices, including the 
model of orbital motion as well as an introduction to the fundamental theory of 
atmospheric modeling. 
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2. Theoretical Background 


This chapter presents the basic equations relating the orbital decay rate of satellites 
to atmospheric density followed by an introduction to the Kalman filter equations. 

Measurement of Drag on Satellite Orbits 

The development of atmospheric models has relied historically on the measurement 
of drag as indicated by the change in period, or the orbital decay parameter (ODP) of 
satellites, and more recently on space-borne instruments such as mass spectrometers and 
gas analyzers. 4 To develop such models though, analytical relationships were required 
that related the orbital decay of a satellite to drag and density. These expressions form the 

basis for the first research approach of this thesis in developing a method to estimate n or 
ODP of LEO satellites by comparing ODP to density predicted by atmospheric models. 

The equations that estimate n are broken down into various categories based on 
satellite orbit type. For simplicity, the atmosphere has been assumed to be spherical and 
exponential, with constant density scale height H, and rotating at the angular rate of the 
Earth. 5 

For satellites with e < 0.2 and ae/H > 3, n is estimated by 6 

n = 3?r(5an : p p exp(--^) I 0 + 2el, +^e 2 (l 0 +I,)+ 0(e 3 ) (1) 

H [ 4 

where p p is the density at perigee, Io, Ii, and L are Bessel functions of the first kind of 
order n and argument Z, written I n (Z), where Z = ae/H. Further, 6 is given by 
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s 


( 2 ) 


FAC d 

m 


where C D is the coefficient of drag, A is the reference area of the satellite, m is the mass of 
the satellite, and F is the atmospheric rotation factor represented by 


F = 


^ r w 
1 cost 

V > 7 


for i < — 


(3) 


F = {2 - F(?r - /)} for i > £ 

2 


where r p and v p are the position and velocity of the satellite at perigee relative to the center 
of the Earth, and w is the angular rate of rotation of the atmosphere. For satellites with 


e < 0.2 and Z < 3, n is approximated by 

n = 3;r<San 2 p p exp(-Z)[l 0 + 2el, + 0(e 2 )] 

However, for satellites in orbits with eccentricities larger than 0.2, n is given by 7 

1 3 

n 1 X „ 2 f Ha V (l + z)~2 ( 8 e- 3 e 2 -l / 0 


(4) 


(5) 


These analytical approximations to satellite n can be used to improve atmospheric models 
if the ballistic coefficient of the satellite is well known, where the ballistic coefficient, (3, is 

defined to be — - , 8 Calculating the ballistic coefficient requires knowing how a satellite 

is oriented throughout its orbit, which enables the cross-sectional area perpendicular to the 
direction of motion to be calculated. For this reason, spherical satellites are often used 
when atmospheric models are being developed. 9 
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When estimating ODP values given an atmospheric model, calculating the ballistic 
coefficient for satellites in the satellite catalog would be a difficult task, due to the many 
complex satellite geometries existing, and the fact that the orientation of the majority of 
satellites is unknown, and perhaps changing. However, given actual ODP values from 
observations and predicted density and perigee, the mean ballistic coefficient of a satellite 
could be estimated by using the appropriate ODP equations. 

Despite the obstacles in evaluating ODP, the feasibility of estimating satellite ODP 
by fitting predicted model densities to actual ODP in a least squares (LS) sense will be 
evaluated in Chapter 5. The degree to which the estimated ODP values agree with actual 
data will depend on the accuracy of the atmospheric model, assuming the particular 
satellite in question is not thrusting or subject to other perturbations other than drag and 
the geopotential. There are several factors that greatly affect density and it is the 
modeling of these factors that will determine how well an atmospheric model will fit the 
ODP of a given satellite. The handling of these factors varies from model to model and 
for this reason, a number of atmospheric models will be applied towards the estimation of 
ODP. But prior to beginning this, a brief discussion of the major factors that affect 
density is introduced in Chapter 3, followed by short descriptions of four atmospheric 
models. 
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Discrete Kalman Filter 


The second research approach, a single state Kalman filter 10 , will be applied to 
update the state or ODP after each fence observation. Kalman filtering is a linear, LS 
process that uses state-space methods and recursive algorithms to estimate a state 
variable, or signal, from measurement data containing an element of random noise. One of 
the main features of Kalman filtering is using the results from the previous step to estimate 
the current state, as opposed to a batch process that generally uses all the data to estimate 
the state. In addition, a batch process can not readily handle state noise, particularly if it is 
time dependent. In the case of estimating ODP of a LEO satellite, it will be shown that 
ODP is highly dependent on the time varying interaction of solar energy with the upper 
atmosphere as indicated by the geomagnetic index, A p , and solar radio flux, F 10 . 7 , which 
will be introduced in Chapter 3. It is the time varying state noise-level of ODP that 
suggests Kalman filtering will produce a better state estimate during periods when the 
state noise-level of ODP is rapidly changing, as is the case during the waxing and waning 
of solar and geomagnetic storms. 

To take advantage of the Kalman filter recursive solution, a few assumptions must 
be made. First, the process to be estimated can be modeled in the following form: 

x k + 1 = i* k + w * ( 6 ) 

where the state vector at time U , denoted by x* , <£>*, *+ 1 , represents the state transition 
matrix (STM), and w k , a white noise sequence with covariance Q* . The STM describes 
how the state evolves from one time to another. White noise is defined as a sequence of 
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random variables uncorrelated in time with a mean of zero. Therefore the covariance 
matrix can be expressed as: 

E [ W * W J] = QA., (7) 

where 5^ is the Kronecker delta function. This means that there is no statistical 
relationship between the value of w* and w* +1 , i.e., for any time other than when t* equals 
t j. However, this does not mean that the individual elements of w* are uncorrelated. 

There could be a significant correlation at any time as described by the off-diagonal terms 
in Q* . The second assumption necessary to utilize discrete Kalman filtering requires the 
observation or measurement of the process to occur at discrete points in time according to 
the linear relationship: 

z k = x* + v t (8) 

Here, z* denotes the vector of measurements at time t*, H* is a matrix defining the ideal 
(noiseless) mapping from state space into measurement space at time t* , and lastly v* 
symbolizes the associated measurement error and similarly to w*, assumed to be a white 
sequence with known covariance structure 

E [ v * v y] = R A.y ( 9 ) 

and having zero cross-correlation with w* that is 

E [ w * v ; ] = 0 ( 10 ) 

for all k and j. A third requirement is that there is information available regarding the 
initial estimate of the state at time t*, and that this information is based on knowledge prior 
to t*. This prior or a priori estimate is represented by , where the symbol A or “hat” 
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above the state vector denotes that this is an estimate, and the minus means that this is a 
best estimate prior to processing the information in measurement z* at time t*. A further 
assumption is that the associated error covariance of is known. Estimation error is 
defined as the difference between the true and estimated state, i.e., 

e * = (id 

The error covariance matrix is then: 

P; = E[e;e; T ] (12) 

where it has been assumed the estimation error has a mean of zero, i.e„ the estimate is 
unbiased. 

With the above assumptions, it is now possible to utilize the measurement z k to 
improve the prior estimate. This is accomplished by" 

xj = x; + K^Zj-H^xJ) ( 13 ) 

where K* is a linear weighting factor applied to the difference between the measurement z k 
and the resultant linear mapping of the prior estimate xj from state space into 
measurement space. The gain K* is determined by picking an optimization criteria. The 
Kalman gain is determined by picking the gain that minimizes the terms along the diagonal 
in the error covariance matrix P* This is chosen because the diagonal terms represent the 
estimation error variances for the state vector elements. The Kalman gain K* is then given 
by: 

K * = P;Hj(H t P;Hj + R, )-' (14) 

The covariance matrix associated with the optimal estimate shown in (13) is then given by 
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(i- K * h* )P; 


(15) 


p; = 

where I is the identity matrix. All the elements are now present to update the estimate and 
the associated error covariance matrix, but to be able to continue the process in 
anticipation of the next observation at l k+l , both x* and P* need to be projected ahead in 
time to function as the next a priori associated with measurement z* +] . This is 
accomplished via the STM. 

*1 (16) 

p; +1 = + q* (i7) 

This process is then repeated as future measurements become available. Figure 2 
illustrates the Kalman filter process. 12 
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Figure 2: Kalman filter process 
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3. Atmospheric Models 


Empirical atmospheric models have been developed for a number of reasons, some 
of which include improving orbit prediction capabilities, understanding atmospheric 
processes and, more recently, aerobraking maneuvers. 13 There are three basic types of 
empirical atmospheric models that have been developed: those developed from studying 
the rate of change in the period of satellites due to drag, those formulated from space- 
borne instruments such as mass spectrometers, and those models which are hybrids or 
comprised of both drag-derived and instrument based measurements. 

To provide background prior to the development of a density based ODP model, 
the input requirements of atmospheric models will be discussed, followed by descriptions 
of the four models of the upper atmosphere that were chosen for this study. An 
introduction to atmospheric modeling can be found in Appendix A. 

Input Requirements 

Table 1 shows the inputs required by the atmospheric models before density can be 
estimated. 
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Table 1 : Required atmospheric model inputs 


| Argument 

T Units/Form "1 

Geodetic Altitude 

km 

Geodetic Latitude 

Degrees 

Longitude 

Degrees 

Day of Month 

DD 

Month 

MM ~ 

Year 

YY 

Hour of LST 

HH 

Minute of LST 

mm 

Previous Day Solar 
Radio Flux. Fio? 

x!0’“ Watts/m 2 /Hz 

81 -Day Smoothed 
Solar Flux, F ]0 7 

xlO" Watts/nr/Hz 

Current Day 
Geomagnetic Index, 



0-400 


Geodetic Altitude 

The density of the atmosphere is primarily a function of the altitude above the 
surface of the Earth, which is oblate. As a result, the atmosphere is also oblate meanin* 
the atmosphere has an equatorial bulge as well.' 4 Therefore, to model the atmosphere 
correctly, a reference must be chosen to represent the surface of the Earth from which 
altitude will be measured. Various models of the surface of the Earth exist, but the most 
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common in use is a reference spheroid, which is an ellipse rotated about the minor axis to 
represent the oblateness of the Earth. Specifically, the ellipse is defined by the equatorial 
radius, R ffi = 6378.140 km , and the ellipticity or flattening. 


f 



1 

298.257 


(18) 


where R p is the polar radius. A much more complex model is the equipotential surface of 
the gravitational field, which is known as the geoia or mean sea level, and contains many 
local irregularities due to the non-uniform mass distribution of the Earth. 15 True geodetic 
altitude is the altitude measured from the geoid upward, but the basic reference spheroid 
model is used to approximate geodetic altitude and is adequate for use in most 
atmospheric models. 

Geopotential Altitude 

Gravitational potential is defined as 

Z 

O = jgdZ (19) 

0 

where <J> is the gravitational potential, g the acceleration due to gravity, and z the 
geometric height above a reference spheroid. Geopotential height or altitude is expressed 
as 
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h 


d> 

8 o 


( 20 ) 


■ “ -\s*z 

So 8 oo 


where h is the geopotential altitude and g 0 is the fixed reference value of gravity, equal to 
9.80665 m/s 2 . 16 


Geodetic Latitude 

Although the oblateness of the Earth does not affect the definition of longitude, it 

does complicate the definition of latitude. Figure 1 illustrates two common definitions of 
latitude. 



Figure 1 : Geocentric and geodetic latitude 
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The angle labeled X' is called geocentric latitude and is defined as the angle between the 


equatorial plane and the position vector R from the geocenter. Geodetic latitude is 
represented by the angle X and is defined as the angle between the equatorial plane and the 
normal to the surface of the reference spheroid. Geodetic latitude is the foundation for 
most maps and, in this case, geodetic latitude would be synonymous with geographic 
latitude. 17 

Solar Radio Flux 

Ultraviolet (UV) solar radiation heats the upper atmosphere through adsorption 
and consists of two components, one associated with the 27-day solar rotation and 
sunspots, and the other related to the 1 1-year solar cycle. Due to energy adsorption, UV 
solar radiation is difficult to measure directly at the surface. However, the extreme 
ultraviolet (EUV) radiation received from the sun has been highly correlated to the 
surrogate index, Fi 0 . 7 , which is a ground-based measurement representing the solar radio 
flux at a wavelength of 10.7 centimeters. 18 The solar radio flux also consists of the short- 
term solar rotation component and the long-term 1 1-year solar cycle. Both of these 
components affect the upper atmosphere differently and must be treated separately. 
Although separate values of these two components of the solar radio flux are not easily 
available, a relationship is used relating the 1 1-year solar cycle to the flux averaged or 
smoothed solar radio flux, Fi 0 7 . 19 This allows the smoothed solar radio flux to be used to 
represent long-term solar cycle effects on the atmosphere. 
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Geomagnetic Index 


Geomagnetic activity is the result of the interaction of the solar wind with the 
Earth s magnetic field. Solar wind kinetic energy is partially adsorbed by the 
magnetosphere, transformed, and eventually dissipated in the magnetic polar regions of 
the atmosphere in the form of heat and ionization. Geomagnetic activity is monitored by 
means of planetary indices such as A p , which like the solar radio flux F l0 7 , is an indirect 
measure of source strength that has a reasonable correlation to observed energy 
dissipation effects. Specifically, A p is a ground-based measurement taken using 
magnetometers placed at various stations . 20 


Model Descriptions 

Four models were studied and their ability to correlate with satellite drag was 
compared. The models examined were the Marshall Engineering Thermosphere (MET) 21 , 
the L.G. Jacchia 1971 (J71) 22 model, the J77 23 model, and lastly, the Mass Spectrometer 
Incoherent Scatter (MSIS) 24 1986 model. A description of each of these models follows. 

Jacchia 71 Model 

The J7 1 model is a revised version of the drag-derived J70 model. 25 Minor 
modifications have been made to numerical coefficients, as well as to the height of the 
homopause. The homopause is the transition region where the model switches from 
hydrostatic equilibrium, governed by the barometric equation, to diffusive equilibrium as 
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represented in the diffusion equation. Changes were based on mass spectrometer and 
EUV absorbtion data at 150 km, suggesting that the concentrations of N 2 and O 2 needed 
to be decreased by 1 6 and 36 percent respectfully, whereas atomic oxygen needed to be 
increased by 37 percent. In addition, the correction made in J70 to the exospheric 
temperature, T» , due to the semiannual variation (S AV) was replaced in J7 1 with a more 
sophisticated correction made directly to density. The J7 1 model dissociates SAW from 
temperature entirely, clearing up many puzzling results in the helium- hydrogen region and 
eliminating the need for adding ad hoc variations for these constituents. 22 

Given the above required model inputs, the process of determining density begins 
by calculating T„. J7 1 accomplishes this by using empirical relations that correct for 
variations due to solar activity, diurnal or day/night effects, latitudinal/seasonal variations, 
and geomagnetic activity. Once T„ is calculated, the result is then applied to Jacchia’s 
empirical temperature profiles, and, together with an expression for the atmospheric mean 
molecular weight, numerical integration is performed on the barometric equation starting 
with a boundary condition at 90 km and integrating up to an altitude of 100 km Above 
this point, the diffusion equation is integrated separately for each atmospheric constituent 
and the partial densities are combined to provide the total density. Once the integration up 
to altitude Z has occurred, corrections are added to the total density for the seasonal 
latitudinal variation (SLV), SAV, and the winter helium “bulge”. 26 

Marshall Engineering Thermospheric Model 

The Marshall engineering thermospheric (MET) model is a drag-derived model 
that is based on J70 but with the SLV and helium bulge corrections taken from the J71 
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model 21 . Other differences between MET and J71 include modifications made to 
coefficients in the temperature profile and the sixth order polynomial used to estimate 
mean molecular weight between 90 and 100 km. In addition, the height of the homopause 
in MET is set to 105 km, whereas in J71 it is lowered to 100 km. 

MSIS-86 Model 

The mass spectrometer incoherent scatter (MSIS) model was developed by A. E. 
Hedin using spacecraft borne mass spectrometers as well as ground-based incoherent 
scatter radar data. This model uses a Bates temperature profile 27 for the upper 
thermosphere and an inverse polynomial for the lower thermosphere, which are both 
functions of geopotential height rather than geodetic height. Substituting geopotential 
height for geodetic height allows exact integration of the barometric equation, assuming a 
constant mean molecular weight and using a boundary condition at 120 km. Exospheric 
temperature as well as other primary quantities are expressed as functions of geographic 

and solar/magnetic parameters using spherical harmonics in latitude and longitude where 
relevant. 24 

Jacchia 77 Model 

The Jacchia 77 model is a major revision of Jacchia’s earlier models. Here he 
incorporates mstrument based mass spectrometric and extreme ultraviolet (EUV) data in 
an effort to improve the representation of individual atmospheric constituents, while using 
satellite drag data to indicate total density. One of Jacchia’s major changes in J77 is in his 
formulation of the effect of geomagnetic activity on temperature and density. In his prior 
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models, corrections were made only to T„ prior to integration. However, J77 
incorporates both an additive term to T„ and a perturbation to the temperature profile, 
necessitating further integration. Other terms are added as well, including a term 
representing the effect of the magnetic disturbance on the height of the homopause and a 
term modeling an equatorial wave that convects from the geomagnetic pole regions 
towards the equator. 23 
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4. Data Description 


The Catalog of Earth satellites is a database containing information on every 
known object greater than or equal to 4 inches in diameter that is currently in orbit. 28 
Objects include active and inactive satellites, rocket bodies (RB’s), and general space 
debris. Currently, the catalog contains over 8,000 satellites, with each catalog record 
containing: a satellite identification number, orbital elements, an epoch or time associated 
with the elements, and an orbital decay parameter (ODP). 29 

One purpose of the catalog is to enable the detection of launches of new satellites, 
and to determine the country of origin, for security reasons. In addition, the US holds a 
treaty agreement with other nations, requiring the US to predict the time and location for 
reentry of space debris, and to notify the appropriate country where landfall is expected. 30 
Further, with an increasing number of satellites and debris comes the increased need for 
collision avoidance. By tracking satellites and updating the catalog, the catalog may be 
used to accomplish all of these tasks. However, it is during periods of the extreme solar 
and geomagnetic storms mentioned previously in Chapter 3, that maintaining the satellite 
catalog becomes more difficult, and therefore the subject of this study. 

In both the ODP predictor model and the Kalman filter approaches to improve the 
maintenance of the satellite catalog, the primary data source will be ODP values listed in 
the catalog for a dozen LEO satellites. But before attempting to implement these 
approaches, an understanding of the data and its source is required. To accomplish this, a 
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description of the radar fence and fence observations is provided, followed by a discussion 
on the calculation of ODP. 

Radar Fence 

The Naval Space Command (NSC) operates a radar fence across the southern 
United States for purposes of helping to maintain the catalog of Earth orbiting objects. 31 
The fence is comprised of 3 transmitters and 6 receiver stations located on an approximate 
great circle with an inclination of 33.6° relative to the equator. Figure 2 shows a map of 
the U.S. where the symbols indicate the location of the six receiver stations. 



Figure 2: NSC radar fence receiver locations. 
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Transmitters 


The three transmitters are located at Gila River, NM, Lake Kickapoo, TX, and 
Jordan Lake, A 1. Continuous wave (CW) illumination is used to provide the maximum 
average power. Under high power conditions, radar coverage is a thin vertical fan of 
width 0.02°. The resolution of the fence is an important factor when considering the 
growmg size of the satellite catalog. To distinguish the reflected signal of one satellite 
from another, there must be an adequate amount of time between signals. However given 
the above fence characteristics, the radar fence would be approximately 350 meters wide 
at an altitude of 1000 km, providing a temporal resolution of 50 ms. 32 With a catalog of 
8,000 satellites, the number of fence crossings is approximately 20,000 per day, or the 
equivalent of a fence-crossing every 4 seconds. Thus, resolving fence-crossing signals 
using a radar fence with a temporal resolution of 50 ms is quite feasible, assuming drag 
effects are insignificant. However, during periods of severe solar storms, changes in drag 
can alter the crossing time of a satellite by more than 10 seconds. It is during these 
perturbed periods that the resolution of fence-crossing signals and subsequent 
identification of satellites can become a problem, and hence the need and interest in finding 
methods to improve the maintenance of the satellite catalog. 

Receivers 

Up to six stations can receive the reflected signal from a satellite passing through 
the fence, when the elevation angle of the satellite is above the local horizon. 

Fundamental measurements are elevation angle or zenith angle, azimuth angle, and time of 
crossing. Radio interferometry is used to measure the angle an incoming signal has with 
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the local vertical and the “time of crossing” is determined by the peak strength of the 
reflected signal. Typical crossing time can be resolved with errors of no more than several 
milliseconds. Position can be determined using two receiver stations and triangulation, 
with an accuracy of 400 meters. 33 

Observations 

Station observations are recorded in two forms, depending on the source of the 
data. NSC records elevation and azimuth angles, whereas the US Space Command 
(USPACOM) records x, y, z position. Typical observations are recorded in the following 
format: satellite crossing time in the form of year, month, day, hour, minute, and second, 
the source of the data, i.e. NSC fence, or USPACOM, the receiver station number, the 
satellite 5-digit ID, and for NSC observations, elevation or zenith angle and azimuth angle. 
Elevation angle is the angle between the local horizontal and the satellite vector, whereas 
zenith angle is the complimentary angle. Azimuth is the angle formed in the horizontal 
plane between the local East vector and the horizontal component of the satellite vector. 

If the source of the observation is USPACOM, the aforementioned angles are replaced 
with x, y, and z position. 34 

Orbital Decay Parameter (ODP) 

The orbital decay parameter (ODP) is produced using several days of fence 
observations in a differential correction (DC) process that determines the updates to the 
seven-element model. Because the number of observations is typically greater than the 
number of elements, no unique solution exists. Therefore the best solution in a least- 
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squares (LS) sense is sought, or equivalently, a solution that causes the sum of the squares 
of the fence crossing residuals to be a minimum. 35 The orbit type of each satellite 
determines actual data spans. For example, a satellite that is experiencing rapid orbital 
decay will require a short fit span of about a day of observations, whereas for a satellite in 
a slowly decaying orbit, fit spans of 3-7 days or longer are used. After recording the 
required amount of data, the DC process yields mean elements with the epoch occurring 
on the last observation time. 36 This process has several consequences. First, the longer 
the fit span, the smoother the resulting ODP will be because this is an average solution 
that is the best fit over the entire span of data. Second, because the time-tag of the 
elements occurs at the time of the last observation, a phase lag will be introduced into the 
ODP of the satellite of about Vi the fit-span. This means that variation in ODP resulting 
from changes in the atmospheric density will appear in element set data with lower 

resolution and occulting from 1.5 to 3.5 days after the actual atmospheric disturbance 
occurred. 
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5. Empirical Fit of Atmospheric Models to ODP 

Ignoring new satellites, collisions between existing ones, and propulsive maneuvers, the 
major cause of uncorrelated targets (UCT's) is due to changes in the density of the atmosphere 
that are not accounted for in the seven element orbit propagation model used to maintain the 
satellite catalog as described in Appendix B. Variations in solar and geomagnetic activity are 
the primary sources for these disturbances in density. The resulting changes in drag can affect 
the along-track orbital position of a LEO satellite anywhere from several rasters to hundreds of 
kilometers per day, depending on the mass-to-area ratio of the satellite and orbit. 37 Along- 
track deviations change the orbital period of a satellite and can be expressed as changes in 
mean motion or ODP. In addition, the ODP derived from DC of elements is a smoothed or 
averaged estimate of satellite drag which will not contain the resolution necessary to be an 
accurate indication of satellite drag during periods of severe atmospheric disturbance. 
Therefore, applying an averaged ODP that also contains a phase lag can lead to large fence 
crossing errors and UCT’s during times of high solar and geomagnetic activity. However, 
including atmospheric modeling as part of the orbit propagation model by developing and 
applying a model of ODP, these atmospheric disturbances and corresponding 
perturbations in density and drag can be accounted for, and potentially lower the number 
of UCT’s. 

The evaluation of ODP is a complex problem, due to the many factors influencing 
the orbital decay of a LEO satellite. This is revealed in Chapter 2 by the orbital decay 
equations (1), (4), and (5), which relate the ODP of a satellite to density at perigee, the 
ballistic coefficient, density scale-height, eccentricity, mean motion, and semi-major axis. 
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The ballistic coefficient, (3, represented by the quantity ^L. is often unknown or varying 

for many satellites in the Earth catalog. For satellites other than spheres, the calculation of 
the cross-sectional area A may not be straightforward, as in the case of a tumbling 
satellite. In addition, the coefficient of drag is a function of satellite shape, altitude, 

surface characteristics, atmospheric composition, speed ratio of incoming particles, and 
solar activity. 38 

Despite these complexities, an ODP model is developed to estimate ballistic 
coefficients and times of strong atmospheric perturbations, which potentially could be 
used to reduce UCT’s. This will be accomplished by applying different atmospheric 
models to estimate density at perigee for a number of satellites over a period of time. The 
density values will then be correlated to satellite ODP and a least squares (LS) process will 
be used to develop a linear relation between density at perigee and ODP. This approach 
was chosen rather than including density as part of the orbit propagation algorithm 
because adding perturbations due to density and drag would require replacing the 
propagation process with an integration technique such as Cowell’s method. 39 

Prior to beginning the correlation of model densities to satellite ODP values, the 
criteria used in selecting satellites for this study is outlined, followed by a discussion on 
the sources of density variation. 

Satellite Selection Criteria 

To test the approach of correlating model densities to satellite ODP, satellites were 
selected that were in non-circular orbits, having a perigee altitude below 1000 km, and 
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daily updates to the element sets. Further, a time period of interest was selected that 
involved significant atmospheric disturbances as indicated by both the solar radio flux Fi 0 7 
and the geomagnetic activity index Ap. 

The period of interest that was chosen to examine the capability of the various 
atmospheric models was the year 1989. During this period, there were three severe 
geomagnetic storms, the first of which was the 3rd largest in the last sixty years. 40 

Upon implementing the above selection criteria. Cosmos 1220 was chosen. Cosmos 
1220 was launched by the USSR on November 4, 1980 and placed into an orbit with an 
eccentricity of 0.02, inclination of 65°, and perigee height of 575 km. 

Cosmos 1220 is believed to be a cylinder of unknown size and originally capable of 
maneuvering. 41 However, since the period of study of this investigation is over 8 years after 
launch, it is probable that Cosmos 1220 was no longer capable of thrusting. Figure 3 shows 
how ODP of Cosmos 1220 varied throughout 1989. The plot reveals many distinct spikes 
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Figure 3: ODP of Cosmos 1220 during 1989. 


in ODP most likely due to changes in atmospheric density. To verify whether ODP 

variation is driven by changes in density, the factors that influence density and hence drag 
must be examined. 
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Sources of Density Variation 


Prior to correlating the ODP of a satellite to predicted density at perigee, it will be 
useful to understand the major sources of density variation. Atmospheric density has 
many sources of variability that greatly effect the drag of a LEO satellite as it passes 
through perigee or the point of closest approach to the Earth. Some of these fluctuations 
are due to variations in altitude, solar activity, solar rotation, geomagnetic activity, diurnal 
or day/night variations, the seasonal-latitudinal variation (SLV) and the semi-annual 
variation. How well a model predicts each of these variations in density is measured by 
the amount of correlation between the density predicted by the model and the ODP of a 
satellite, and therefore will be used to rank the quality of each atmospheric model. 

Solar and Geomagnetic Activity 

A major source of variation in upper atmospheric density is caused by fluctuations 
of the solar EUV flux and solar wind. To help understand how these factors affect 
satellite drag, the daily solar 10.7 cm radio flux F; 0 7 , and the daily geomagnetic index A p , 
were plotted for the year 1989 in Figure 4. 
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Figure 4. Daily solar flux and geomagnetic activity during the year 1989 


The upper plot displays geomagnetic activity revealing several distinct spikes in A p on a 
background of lower amplitude activity. The spikes illustrate the three significant 
geomagnetic storms that occurred during 1989. The lower plot of the solar radio flux 
shows quasi-periodic” variations that correspond to the sun’s 27-day solar rotation. 
These trends in the solar and geomagnetic activity should be clearly visible in the ODP 
data and also the density data output from an atmospheric model if that model contains an 
adequate representation for such phenomenon. 
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To test this assumption. Figure 5 was constructed as a combination plot of 
atmospheric density indicated by the MET 21 model, ODP of Cosmos 1220, geomagnetic 
index A p , and solar radio flux F | 0 .7 during 1989. 




0 1 1 1 1 1 1 1 1 
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Day of Year, 1989 

Figure 5: MET density, ODP of Cosmos 1220, Ap, and F 10 . 7 . 


The various trends of both the geomagnetic effect and solar activity can be clearly seen in 
both the ODP of Cosmos 1 220 and the density predicted by MET. For instance, in the 
early part of 1989 there were a series of three broad peaks in the solar radio flux, reaching 
well over 250 xlO' 22 W m' 2 Hz' 1 . The third peak was accompanied by the third largest 
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storm in the past sixty years, occurring on the 72nd day of the year. 40 The plot of 
Cosmos 1220 ODP and predicted density show three corresponding peaks during this 
period. Specific spikes and trends in the ODP curve can be attributed to the appropriate 
factor by directly comparing plots of those factors with the plot of ODP. 

Diurnal, Seasonal-Latitudinal, and Altitude Effects 

To help illustrate diurnal or day/night variational effects in satellite ODP plots, 
local solar time (LST) and latitude of perigee were calculated. Figure 6 shows plots of 
LST and latitude of perigee for Cosmos 1220. 




34 





The top plot reveals that the LST' of perigee for Cosmos 1220 goes through several 

diurnal cycles throughout the year. This is mainly due to the precession of the orbital 

plane caused by the oblateness of the Earth. The latitude of perigee plot reveals that 

perigee for Cosmos 1220 is changing slowly. This is because the inclination of Cosmos 

• • 

1220 is near the critical inclination of 63.4° where co approaches zero. With a small co, 
the resultant latitude of perigee will vary slowly. 42 

To assist in interpreting the effect that altitude variation had on density and ODP 
for Cosmos 1220, the height of perigee and the density scale height, H, were calculated. 




Figure 7: Perigee altitude and density scale height for Cosmos 1220. 


f LST equals zero at midnight. 
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with results shown in Figure 7. The top plot shows that the perigee altitude of Cosmos 
1220 increases gradually, with a variation of 20 km throughout the year. Variation in 
perigee altitude occurs primarily due to the oblateness of the Earth. The bottom plot 
displays the values for the density scale height, and reveals that H is fairly constant, with 
an average of approximately 85 km. Because the change in perigee altitude for Cosmos 
1 220 ts only about 25% of the density scale height, altitude changes will represent 
approximately 22% of the total changes in ODP for Cosmos 1220, meaning that the 

majority of the changes in density and ODP will be the result of variations in other factors 
such as solar flux and geomagnetic activity. 


Correlation of Models to ODP 


To begin the comparison of atmospheric models to Cosmos 1220 ODP, model 
inputs were calculated over a period of a year from daily satellite orbital elements. The 
National Oceanic and Atmospheric Administration (NOAA) supplied solar and 
geomagnetic inputs. Atmospheric density was then calculated at perigee over the entire 
year for each of the models described in Chapter 3. Once this was completed, a LS fit 

was performed fitting daily, predicted density to the daily change in mean motion using the 
following linear form: 


ODP,. = A p, (21 

This form was chosen for its simplicity, and due to ODP, i.e., n , being linearly related to p as 

shown in the orbital decay equations (1), (4), and (5). In addition, ODP should approach 
zero when density approaches zero. 
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MET Model 


Figure 8 shows a scatter plot of Cosmos 1220 ODP relative to the linear estimate 
calculated using MET densities in (21). The plot reveals that the relationship between ODP 
and density is mostly linear, however, there is considerable noise present. 
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Figure 8: LS fit of MET density to the ODP of Cosmos 1220. 

The upper part of Figure 9 shows a plot of density as calculated by the MET model using 
Cosmos 1220 element sets. The lower plot shows Cosmos 1220 ODP from element sets, 
(solid line), compared to the linear approximation to ODP using (21) with MET density at 
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perigee (dash-dot). Clearly there is good correlation, but a phase-lag exists between the 
actual ODP and that predicted by MET as indicated by the ODP spikes of Cosmos 1220 



Figure 9: MET density fit showing a phase lag with ODP of Cosmos 1220. 

appearing to the right of those estimated by MET. The plot of ODP also indicates a 
process bias exists and/or an element of noise is present, as suggested by the spikes in the 
negative direction that typically follow large corrections to ODP made by the DC process. 
This may be due to over-correction of ODP. In addition, the increases in density 
associated with the two geomagnetic storms that occurred on DOY 72 and 293 are fairly 
well predicted, yet the density effects from the storm that occurred on DOY 322 are 
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extremely under-predicted. This reveals the geomagnetic effect of the MET model most 
likely needs additional development. But despite the lag, noise, and under-estimation, 
MET generated densities have a 0.64 correlation with ODP from the element sets. 
Recalling the method that was used to generate the ODP contained in the element sets 
discussed in Chapter 4, the existence of a lag between the atmospheric disturbances in 
density and ODP should be expected. 

To remove this lag and reduce the process bias contained in the ODP values, 
running averages were calculated for ODP and density using consecutive data points and a 
phase-shift was introduced into the LS expression in the following form 

ODP, = Cp {l _ k) (22) 

where represents the previous density value at perigee, approximately £-days earlier 
than p,. The actual amount of phase shift was dependent upon the frequency of the DC 
updates to the elements. For Cosmos 1220, updates to the elements occurred an average 
of 1.01 days apart. Thus a fixed phase shift in time of an integer number of days was not 
possible, and hence the need to express this temporal shift as a k - day phase shift, where 
the value k of a satellite is approximated by the average time between updates to the 
elements. Therefore, for the case of Cosmos 1220, k would be approximately 1.01 days. 
Introducing the phase shift in density will counter the apparent temporal lag in ODP and 
should result in a higher correlation. Ideally, density should have been averaged over the 
entire data fit-span to match the smoothed ODP from the differential correction (DC) 
process, but because the fit-spans were unknown, averaging density over the data span 
was not possible. 
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Figure 10 shows the results for Cosmos 1220 after performing a LS fit using a 
running mean of both density and ODP along with a zero and a ljfc-day phase-shift in 
density as calculated by the MET model. The top plot of Figure 10 shows the effect on 
the correlation using running means of consecutive ODP and density values. The bottom 

plot shows the additional increase in correlation due to the introduction of a U-day phase 
shift in density. 



There is still a lag visible between the density and ODP. Therefore, the lag in density was 
increased to 2k and 3^-days with the results shown in Figure 11. The top plot shows the 
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results corresponding to a phase-lag of 2&-days whereas the bottom reveals the effect of a 
phase lag of 3/c-days. 



x 10 -4 



Since the correlation decreased for a phase shift of 3k-days, the above plots suggest that 
the optimum phase shift in density for Cosmos 1220 lies between 1 k and 2/t-days. This 
can be illustrated by plotting the variation in correlation versus the amount of lag applied 
to the density model as shown in Figure 12. 
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Figure 12: MET correlation vs. density lag. 


The existence of a lag in ODP will have a serious impact on fence crossing 
predictions most notably during a period of severe atmospheric disturbance. This 
disturbance in density and drag will not be directly input into the orbit prediction model 
because the orbit model does not contain an atmospheric model. The only way this 
change in density and drag can manifest itself into the orbit predictor model is through DC 
of ODP. However, recall the discussion in Chapter 4 regarding the processing of lengthy 
data-spans of radar fence observations. The updates to ODP from the DC process will be 
an averaged solution, that is the “best fit” over the entire span of observations and lagging 
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up to several days, depending on the fit-span length of the data. Therefore, the DC 
process can not provide the sharp corrections necessary to ODP during periods of severe 
atmospheric disturbance. The end result will be large fence crossing errors during extreme 
solar storms, and thus UCT’s. 

However, applying the ODP model just developed without using a phase lag 
would eliminate the lag between the orbit model and the atmospheric perturbations. The 
phase lag is a process bias formed as a result of applying the DC method to ODP over a 
length of data, and was introduced into the ODP model to match this process bias and get 
a more representative correlation between model densities and observed ODP values. By 
applying the ODP model using atmospheric density predicted at the time of perigee, a lag 
will not exist between density and the estimate for ODP, which then could be used in 
conjunction with the existing orbit model to potentially reduce the number of UCT’s 
occurring during severe solar storms. 

MSIS-86 Model 

Similar results were found for the MSIS-86 atmospheric model, but differences in 
model structure are apparent. Results are shown in Figure 13 and Figure 14 where the 
phase shifting techniques were applied. As previously, the solid line is Cosmos 1220 ODP, 
whereas the dash-dot line represents the estimated ODP calculated using a LS fit of the 
density at perigee predicted by MSIS-86 to actual Cosmos 1220 ODP. 
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The results are similar as with the MET model, but there are several time periods where 
the MSIS-86 model estimates for ODP are lower than MET estimates, as occurs on DOY 
170 and 260. On the contrary, MSIS-86 predicts higher density than MET for the 
geomagnetic storm on DOY 293, suggesting that the MSIS-86 model is representing the 
geomagnetic effect on density more accurately than MET. However, despite the periods 
of greater density prediction, it is apparent that density under prediction is more prevalent, 
given the lower correlation of MSIS-86 to Cosmos 1220. 
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Under or over-estimation of ODP is undesirable. An appropriate atmospheric 
model must be able to properly estimate the major contributors to density variation as 
outlined earlier in this chapter. A model that has diurnal and seasonal latitudinal variations 
represented, but does not estimate the solar flux or geomagnetic impacts to density 
accurately, will not be a very useful atmospheric model for the purposes of eliminating 
UCT’s caused by solar storms. 



Day of Year, 1989 

Figure 14: Correlation of MSIS-86 to Cosmos 1220 using phase shifts of 2 k and 3£-days. 
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Referring back to Figure 4 which shows the solar flux and geomagnetic activity 
during 1989, the time periods of under-estimated MSIS-86 density coincide with periods 
of high solar activity, suggesting that the coefficients governing the solar flux 
perturbations to MSIS-86 density need adjusting. 


J71 Model 

Figure 15 and Figure 16 show the results of the LS fit of the J71 22 model to 
Cosmos 1 220 ODP for various phase shifts in density. 


x 10~ 4 
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Very similar results to the MET model can be seen. However, notice that the J7 1 estimate 


for ODP associated with the A p storm on DOY 293 is, as in the MSIS-86 case, greater 
than the corresponding MET estimate for ODP, implying that the disturbances in density 
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Figure 16: Correlation of J71 to Cosmos 1220 using phase shifts of 2k and 3&-days. 

that occurred on this day due to the geomagnetic effect may be more accurately 
represented by the J71 and MSIS-86 models. 
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J77 Model 


The Jacchia 77 model was evaluated using the same procedure as for the other 
models, with results shown in Figure 17 and Figure 18. 



A general inspection of the J77 estimate for Cosmos 1220 ODP suggests that J77 does not 
predict total density as well as the other models. The spikes in ODP due to the three 
severe geomagnetic storms that occurred on DOY 72, 293, and 322 are all under- 
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estimated, considerably more than any of the other models. 

Analysis of Model Performance 

In comparing the responses of all four models, differences arise in how well they 
match fluctuations in ODP. The majority of these disturbances in ODP are most likely 
caused by variations in the solar flux and geomagnetic activity, as suggested earlier in the 
sub-section entitled “Geomagnetic and Solar Activity”. Comparing Figure 1 1 and Figure 
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14 showing the correlation of the MET and MSIS-86 models to ODP respectively, reveals 
that the MET model responds stronger and with better resolution to variations in solar 
flux than does the MSIS-86 model in this particular case. Another interesting region lies 
in between days 170 and 225 of the year. Both the beginning and the end of this period are 
characterized by strong increases in the daily solar flux. The result forms a semi-circular 
shape in the ODP curve. The MET model responds strongly to the F 10 . 7 variation both at 
the 160-day mark and the 225-day mark. On the other hand, the MSIS-86 model shows 
very little response at the 170-day mark, but estimates the 225-day peak fairly well. 

The response of J71 to geomagnetic activity is most similar to MSIS-86 
predictions. However, J71 estimates for density are slightly higher as can be seen by 
comparing Figure 14 and Figure 16 on DOY 72 and 293. For changes in density due to 
solar flux changes, J71 and MET are most alike. Nonetheless, in comparing Figure 1 1 and 
Figure 16, it can be seen that J71 predicts lower density than MET on days 15, 40, and 
170. But on day 225, the difference is only slight, and on days 125 & 290, J71 is 
considerably larger than MET. There is an apparent difference between the two models 
that seems to appear semi-annually. This phenomenon behaves similarly to the semi- 
annual variation (SAV) which generally peaks twice a year, once in the spring and once in 
the fall. SAV implementation differs between these two models and could be the reason 
for the observed contrast. In MET, the change in density due to SAV is represented by a 
correction to exospheric temperature, T», that is a function of time of year and solar radio 
flux. In contrast, J71 represents SAV by a direct correction to density that is a function of 
altitude and time of year, but eliminates any functional dependency on solar flux. 43 
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The models were fair when it came to reacting to geomagnetic events. Recall 
Figure 4 showing the geomagnetic and solar activity throughout 1989. There were three 
main geomagnetic events, occurring on days 72, 293, and 322. All models responded 
fairly well to the first two geomagnetic storms by correctly estimating a sharp rise in ODP. 
However, the last storm was extremely under predicted by all of the models. Key 
variables for determining atmospheric density were identified for the three storms in an 
effort to isolate a possible cause or deficiency in the models. The severity of the three 
storms happened to decrease for days later in the year, meaning the order of the storms in 
decreasing strength was day 72, 293, and 322. Upon inspection of the ODP plot, the 
order appeared to be contradictory, since the largest spike in ODP occurred during the 
weakest storm. In contrast, it would be expected that the largest geomagnetic disturbance 
would create the greater increase in density, all else being equal. However, after 
examining other variables it was determined that perigee for Cosmos 1220 occurred at low 
latitude during the largest geomagnetic event and was located at 65° and 60° for the other 
storms, respectively. This would help explain why ODP was smaller for the first storm 
when Ap was 249, compared to the third storm with an Ap of 138, since Cosmos 1220 
perigee occurred close to the north geomagnetic pole of the Earth for the third storm*. It 
is in the magnetic polar regions of the thermosphere where thermal and density effects 
from geomagnetic activity would be greatest. 44 Nonetheless, this does not explain why 
ODP for the second largest storm on day 293 was smaller than the ODP that took place 
during the weaker storm on day 322, since perigee was even closer to the north 


J The north geomagnetic pole occurs at a geodetic latitude of 78.6° and a longitude of 289.3°. 
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geomagnetic pole during the stronger storm on day 293. If LST of perigee is considered, 
the resultant behavior of ODP is explained. During the weaker storm on day 322, LST 
was approximately 3 AM, whereas during the stronger storm on day 293, the LST was 
approximately 12 Noon. Geomagnetic activity affects atmospheric density the most 
during the night when solar energy input is at a minimum and therefore temperature and 
density reach their minimums as well. Thus any addition of energy during the night is 
going to have a larger impact on temperature and density and therefore ODP, than an 
equivalent storm occurring during the day. 

Having reasonably explained the pattern of ODP for the geomagnetic storms, that 
still leaves understanding why all the models failed to estimate ODP correctly during the 
third magnetic storm. One possibility is the fact that none of the Jacchia models evaluated 
include LST or longitude in the functions that represent the changes in density due to the 
geomagnetic effect. LST and longitude appear to be important variables, given the 
difference between the effects of the geomagnetic storms on ODP during days 293 and 
322 and the LST when they occurred. The J77 model contains the most complex 
geomagnetic modeling of all the Jacchia models, yet failed to outperform the other 
models. However, Hedin does use LST when factoring in geomagnetic effects in MSIS- 
86, but it is possible that some coefficients need adjusting. 

Extended Satellite Study 

To further test the approach of modeling satellite ODP using atmospheric density 
and continue the evaluation of the 4 atmospheric models, the study was expanded to 
include eleven additional satellites in various types of orbits during the year 1991. 
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Elevated geomagnetic activity during this period was generally more frequent than that 
during 1989, although the magnitudes of the storms were less severe. Table 2 summarizes 
the properties of the satellites, if known, including the 5-digit satellite identification 
number, international reference number, name, height of apogee and perigee, inclination, 
radar cross-section (RCS), shape, size, and mass. 


Table 2: Satellite characteristics 


Sat. 

No. 

Int. 

No. 

Name 

H a 

X 

Hp 

(km) 

i 

(deg.) 

RCS 

(m 2 ) 

Shape 

Size 

(m) 

Mass 

(kg) 

60 

1960- 

014A 

Explorer 8 

1498 
x 399 

50 

0.47 


0.76 L x 
0.76 Dia. 

41 

614 

1963- 

025B 

Hitch- 
Hiker 1 

2759 
x 330 

82 


m 

0.30 L x 
0.90 Dia. 

79.8 

1616 

1965- 
07 8 A 

Atlas D 
RB 

m 

144 

1.62 

Cyl. 

2.05 L x 
0.72 Dia. 

70 ? 

2389 

1966- 
70 A 

OV3-3 

3342 
x 354 

81 

1.01 

Oct. 

0.74 L x 
0.74 Dia. 

75 

2404 

1966- 

70B 

OV3-3 RB 

1286 
x 318 



Cyl. 

1.5 Lx 
0.46 Dia. 

24 

3342 

1968- 

066D 

Explorer 
39 Debris 

1761 
x 619 

81 

0.01 

? 

? 

7 

4222 


Scout B 
RB 

1621 
x 364 

103 

1.60 

Cyl. 

1.5 L x 
0.46 Dia. 

24 

8368 

1975- 

100C 

DELTA 1 

RB(2) 

a 

23 

2.85 

Sphere 

-Cone 

1.32 Dia. 
to 0.94 
Dia. 

66 

11791 

■ 

I 

Atlas F 
RB 

12215 

x264 

63 

2.27 

Cone- 

Cyl. 

1.85 L? 
0.63-1.65 
Dia. ? 

163 ? 

12069 

1980- 

087B 

Atlas 

Centaur 

RB 

9336 
x 268 

■ 


Cyl. 

8.6 Lx 
3.0 Dia 

1815 

15679 

■ 

■»«5B 

Ariane 3 
RB(3) 

30156 
x 264 

6 

15.9 

Cyl. 

9.9 Lx 
2.6 Dia. 

2150 

(e) 
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Results of the extended satellite study are shown in Table 3, which lists the satellite 
number, name, average time k between element set updates, the amount density was 
shifted in &-days, and the resulting correlations of the atmospheric models to the ODP data 
of the above eleven satellites. 


Table 3: Model correlation coefficients 


Sat. No. 

Name 

k 

(days) 

Lag 

(A>days) 

MET 

MSIS 

J71 

377 I 

60 

Explorer 8 

0.75 

3 

0.55 

0.62 

0.58 

0.50 

614 

HitchHiker 

1 

0.84 

3 

0.82 

0.86 

0.88 

0.81 

1616 

Atlas D RB 

0.79 

3 

0.72 

0.74 

0.80 

0.74 

2389 

OV3-3 

0.85 

3 

0.75 

0.74 

0.83 

0.78 

2404* 

OV3-3 RB 

0.74 

3 

0.74 

0.77 

0.77 

0.74 

3342 

Explorer 
39 Debris 

2.21 

1 

0.73 

0.68 

0.71 

r 0.79 



0.69 

3 

0.84 

0.86 

0.87 

0.83 

8368 


0.66 

3 

0.62 

0.63 

0.75 

0.64 

11791 

Atlas F RB 

1.87 

2 

0.92 

0.88 

0.93 

0.85 

12069 

Atlas 

Centaur 

RB 

1.58 

2 

6.71 

0.67 

0.82 

0.72 

15679* 

Ariane 3 
RB(3) 

2.11 

2 

0.74 

0.79 

0.72 

6.57 


* Correlation timespan was limited to DOY = 1-250, due to orbital decay. 

* Correlation timespan was limited to DOY = 1-200, due to orbital decay. 
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The results reveal that the J7 1 model performed the best out of the 4 models as 
indicated by the correlation of J71 being the highest for 7 of the 1 1 satellites. This 
suggests that J7 1 is the better atmospheric model for the altitude range studied. In 
general, it can be seen from the correlations that all the models performed well for nearly 
all of the objects, having correlation coefficients from 0.7 to 0.9. The only exception is 
Explorer 8, where the correlations are in the range of 0.5 to 0.6. There is a myriad of 
reasons why the models didn’t fit this satellite as well as the others. Examining a plot of 
the correlation of MET density to ODP of Explorer 8 in Figure 19 reveals that the lag 


x ^ 0~ 4 



Figure 19: Correlation of MET to Explorer 8 using phase shifts of 3&-days. 
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ODP, Revs/Day A 2 


between the ODP data and the density model still exists and is more than likely the reason 
for the low correlation. The presence of a lag even after shifting the density 3k-days, 
where k for Explorer 8 is 0.75, suggests that the data fit-span used in calculating the ODP 
of Explorer 8 is longer than 6*-days, i.e„ 6 x 0.75, or approximately 4.5 days. Investigating 
further, the density was shifted out to a total of lOk-days. Figure 20 shows a plot of the results 


of shifting the density 8k-days. 
x 10 -4 
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As in previous plots, the solid line represents actual ODP, whereas the dash-dot indicates 
the ODP estimated by the MET model. By shifting the density, the lag between the averaged 
ODP solution from DC and the atmospheric model has been eliminated. The end result is 
an increase of 25% in the model correlation from a value of 0.55 up to 0.69. Figure 20 
shows the lag versus correlation for Explorer 8. 



Figure 21: Effect of lag on Explorer 8 correlation 

It appears that the optimum lag falls between 8 to 9&-days, which is equivalent to 8.5 x 
0.75, or a lag of 6.4 days. Recall that the timestamp from the DC process occurs at the 
time of the last observation so that the lag will be approximately V 2 of the fit-span. 
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ODP, Revs/Day A 2 


Therefore, a lag of 6.4 days would suggest a fit-span of about 13 days. Explorer 8 has 
been orbiting for over 35 years due most likely to its relatively high height of perigee and 
small RCS, which translates into lower drag. Hence the orbit of Explorer 8 is decaying 
very slowly. A slowly decaying orbit does not require frequent DC of elements as would 
be required for a rapidly decaying orbit. Therefore, observations of a slowly decaying 
orbit such as Explorer 8, can be extended over many days for a longer fit-span, and thus, 
lower the computational burden associated with performing frequent DC updates. 

Another satellite of interest is the Atlas F rocket body (RB). Correlations for this 
satellite ranged from 0.85 to 0.93. To get a better understanding as to why the 



Day of Year, 1991 

Figure 22: Correlation of J71 to Atlas F using phase shifts of 3k-days. 
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correlation was significantly higher than the other satellites, a plot of the LS fit of J7 1 to 
the ODP of Atlas F was examined as shown in Figure 22. The plot shows that the ODP of 
Atlas F has a large long-term linear trend superimposed upon finer short-term variations 
caused most likely by changes in atmospheric density resulting from fluctuations in solar flux 
and geomagnetic activity. The linear trend in ODP for Atlas F is most likely due to the 40 km 
decrease in the height of perigee. Large decreases in altitude equate to increases in density and 
ODP. When the change in altitude becomes large enough, i.e. approaching the value for the 
density scale height, H, ODP becomes a function dominated by altitude changes rather than by 
variations in solar and geomagnetic inputs, which appears to be the case for Atlas F. All of the 
models were able to match this large linear tendency in ODP, which accounts for the high 
correlation. However, it is more critical that the models match the finer details of the ODP 
curve, which would indicate that the models are following the physical short-term trends 
occurring in the atmosphere. For this reason, the high correlation given by the models for 
Atlas F should be viewed with caution. Ideally, the large linear trend should be removed before 
analyzing the performances of the models. This can be accomplished by passing the ODP of 
Atlas F through a high pass filter. A high pass filter removes low frequency content from 
the input signal while allowing the high frequency signal to pass through unaltered. 45 

As mentioned earlier under satellite selection criteria, it is desirable to study 
atmospheric drag effects using satellites in fairly eccentric orbits, where the drag effect would 
occur primarily near perigee, as opposed to a circular or nearly circular orbit where the altitude 
of the satellite would be nearly the same throughout the orbit, and thus latitudinal, seasonal, 
and LST variation of drag and density could not be isolated. It should be expected that for the 
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circular orbit case, ODP would not be a strong function of H, since altitude is nearly fixed. 

However, for eccentric orbits, altitude varies and therefore variation in H should have a 
stronger effect on ODP. 

Recall the relationship expressed in (5) where satellite ODP, or n , was proportional 

1 

to pH for orbits with eccentricities of 0.2 and larger, and therefore could be used to 
estimate the ODP of a satellite using densities predicted by an atmospheric model. Of the 
eleven satellites examined, a number of them were in orbits of eccentricity greater than 

0.2, and given this relationship, there should be an improvement in model correlation if 

1 

the expression pH 2 is related to ODP instead of simply p. To investigate, pH 2 was 
calculated for all the satellites and correlations performed to their respective ODP values. 

Table 4 shows the percent change in correlation for the MET model using pH^ relative to 
the correlations and phase lags given previously in Table 3. 
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Table 4: Change in correlation after fitting pH 3 to ODP 


| Sat, No. 

Name 

Change in Corn ( %) 

| Eccentricity 

60 

Explorer 8 

-1.60 

0.075 

614 

HitchHiker 1 

-0.60 

0.150 

1616 

Atlas D RB 

+0.14 

0.150 

2389 

OV3-3 

+0.40 

0.180 

2404 

OV3-3 RB 

+4.44 

0.060 

3342 

Explorer 39 
Debris 

-0.40 

0.073 

4222 

Scout B RB 

+0.24 

0.085 

8368 

DELTA 1 RB(2) 

+8.80 

0.330 

11791 

Atlas F RB 

+0.55 

0.470 

12069 

Atlas Centaur RB 

+3.94 

0.400 

15679 

Ariane 3 RB(3) 

+2.17 

0.690 


As was suspected, correlations increased for those satellites in orbits with eccentricities 
larger than 0.2. For most of the satellites with e < 0.2, the correlation coefficients went up 
only slightly or in some cases decreased. There were exceptions however, notably the 
OV3 3rd stage RB, where the orbit was fairly circular yet the correlation increased. This 
was most likely due to decay of the OV3 orbit, causing perigee altitude to decrease into a 
denser atmosphere where H changes more rapidly and therefore strengthening the 
dependency of ODP on H. 
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Ballistic Coefficient Estimation 


Ballistic coefficients were estimated using predicted densities at perigee in 
conjunction with the appropriate theoretical orbital decay equation (1), (4), or (5) from 
Chapter 2. The empirical ballistic coefficient of a given satellite was obtained by 
equating the theoretical ODP expression to actual ODP values, and solving for p. A 
theoretical p was then calculated if possible, based on available satellite physical 
characteristics, i. e . , shape, size, and mass. Because the cross-sectional area perpendicular 
to the velocity direction of non-spherical satellites can vary throughout an orbit, as is the 
case for a tumbling satellite, minimum, maximum, and mean areas were calculated, and the 
results used to estimate the possible ranges in satellite P‘s. 47 

The results of estimating ballistic coefficients for all the satellites are shown in 
Figure 23, where the circles represent P as estimated by MET, and the asterisks are a 
mean P as calculated from satellite physical characteristics published by King-Hele. 29 
Error-bars appear depicting the possible range in P based on maximum and minimum 
satellite reference areas if known, for the non-spherical satellites. 
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Satellite ID 


Figure 23: Ballistic coefficient plot. 

The plot shows the results for the satellites from the extended satellite study and Cosmos 
1220, and reveals the good agreement between the empirical and theoretical estimates for 
p. Notice that for satellite 3342, also known as Explorer 39 debris, (3 is extremely small. 
No physical characteristics for Explorer 39 debris are available, therefore a theoretical P 
could not be calculated. 
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In summary, the atmospheric models correlated fairly well to satellite ODP data, 
with the J71 model performing the best. By applying the LS fitting process between 
model densities and satellite ODP values, two simple linear ODP models were developed 
that could be used to predict satellite orbital decay rates. The first form where ODP was 
estimated by a LS fit to density, would be used for quasi-stable orbits with e < 0.2, 

whereas in the second form, ODP was fit to pH^. This version would be reserved for 
orbits with e > 0.2 and for satellites that are experiencing high rates of orbital decay prior 
to re-entry. In the latter case, the density scale height, H, would be more likely to change 
significantly enough to effect satellite orbital decay rates and hence must be accounted for 
in the ODP model. In addition, it was determined that a phase lag exists between satellite 
ODP values as determined by DC, and atmospheric disturbances. The amount of lag 
depends on the fit-span that is used to perform the DC. How well the atmospheric models 
predicted ODP fluctuations brought on by solar disturbances determined the level of 
correlation. Although the models generally represented the atmospheric dynamics well, 
there were occasional instances where all the models underestimated effects of 
geomagnetic and solar activity, clearly indicating the need for further advancement in 
atmospheric modeling. By using the above ODP models and eliminating the lag, improved 
fence crossing predictions could be made, potentially reducing the number of UCT’s. In 
addition, atmospheric models can be used in conjunction with orbital decav theory to 
determine satellite ballistic coefficients, which can aid in identifying satellite class or be 
used to monitor changes in (3 that might occur due to alterations in the orientation of the 
satellite to the orbit-plane, or changes in the mass-to-area ratio due to maneuvers or 
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collisions, etc. The empirical ballistic coefficients of the majority of satellites agreed well 
compared to their theoretical values. 



6. One-State Kalman Filtering Approach 

Recall from Chapter 4 that ODP from element set data is the orbital decay 
parameter calculated from a DC process of variable data length. This makes the ODP 
smooth and lag the physical processes occurring in the atmosphere. Using smoothed ODP 
data with a lag to make fence crossing predictions can lead to large crossing errors during 
severe geomagnetic storms such as the one that occurred in March of 1989, when crossing 
errors grew to over 10 seconds for certain satellites. To circumvent this, a process must 
be developed that would provide a more responsive ODP with little or no lag. 

The purpose of this chapter is to determine the feasibility of using a one-state 
discrete Kalman filter to process observations, i.e., take in measurements as they occur at 
the fence and make optimal corrections to ODP, i.e., the state, as required to reduce the 
number of UCT’s occurring during large solar storms. Kalman filtering had been applied 
previously for the same purpose, but a three-state filter had been utilized. 48 

A one-state filter was chosen because the recursive filter equations listed in 
Chapter 2 are reduced from vectors and matrices to simpler scalar expressions. The result 
is a filter that is easier to implement and maintain and therefore more readily applied to a 
large-scale operation such as maintaining a satellite catalog, as opposed to using a multiple 
state filter. However, the disadvantage of modeling the physical world by a single state 
variable is that the model will be more likely to be deficient in representing real world 
dynamics. Thus the more state variables there are describing state dynamics, the more 
accurate the model will be, but at the price of adding more complexity to the operation. 48 
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State Model 

The single state variable that was chosen for filtering was ODP, i.e., n , making the 
state “vector”: 

x = [ODP] (23) 

Normally as part of the Kalman filter process, state dynamics are handled by the STM cp, 
where dependencies among state variables are expressed using partial derivatives. As 
stated previously, ODP affects both i and n as seen in equations (50) and (51). Since 
ODP is the only state variable, the state transition matrix (STM) will simply be unity. This 

means that changes in i and n due to changes in n can not be reflected by applying the 
STM to the state. To circumvent this, changes in n and i due to variation in the state 

variable can be reflected by repropagating n and t in time, once the new estimate for n 
or ODP is determined. 

State Noise Model 

State noise is a measure of the uncertainty in the physical variable being modeled. 

In this application, state noise would be the residual variations in ODP as indicated by the 
existence of fence-crossing errors. If the current value of ODP is predicting crossing 
times with zero error, then there would be no uncertainty in the knowledge of the state 
and hence state noise would be zero. However, this is rarely the case since ODP is 
constantly changing primarily due to atmospheric drag for LEO satellites. 

Modeling the state noise associated with a dynamic atmosphere can be a 
challenging task. Recall earlier the strong perturbations to atmospheric density and hence 
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ODP that caused by fluctuations in both the solar flux, F I07 , and geomagnetic index, 

A p Given the dynamic nature of atmosphere, it can be expected that the noise levels of 
the atmosphere are highly dependent on the values of F l0 . 7 and A p . Thus it can be 

expected that state noise levels will be quite different for quiet atmospheric conditions 
where F 10 . 7 < 100 and Ap < 50, as opposed to “noisy” periods during severe solar storms 
where F 10 . 7 > 200 and/or Ap > 100. This suggests that a representative state noise model 
must be able to adapt to the changing conditions of the atmosphere, in contrast to a state 
noise model that uses a fixed noise level based on the average value of ODP in the past. 

To investigate which type of process noise model would work best, both models 
were evaluated using real fence data for the satellites listed in Table 5. First, the fixed 
noise model used was of the form 


Q*= [r dt ] 2 


(24) 


where dt is the time between observations in days and y is some constant to be 
determined. Second, the adaptive noise model used was of the form 

Q* = [a ODP<_, dt] 2 


where a ODP*. , is a percentage of the previous state estimate ODP*. , and is used to 
represent the current level of state noise. In addition, Q* in both expressions is 
proportional to dt 2 . Thus the state noise increases with the time between fence 
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observations, due to the higher uncertainty in older estimates of ODP representing the 
current rate of the orbital decay of a satellite. 

Using ODP directly in the state noise expression allows the state noise model to be 
adaptive to the dynamics of the atmosphere. That is, if a perturbation of some sort occurs 
in the atmosphere, it will manifest itself in the estimated value of ODP output from the 
Kalman filter. This in turn will affect the state noise level, increasing noise for perturbed 
periods whereas decreasing it for quiet intervals. However, there is the possibility that 
measurement errors exist, e.g., associating a fence-crossing observation with the incorrect 
satellite, errors in data transmission can occur, etc. If the low quality of the measurement, 
e.g., at time t k , was not indicated by a corresponding rise in the measurement noise R k , 
and the error covariance P k .| has converged, the Kalman filter will react to the noisy 
measurement by adjusting the state estimate of ODP to eliminate the residual crossing 
error and result in an erroneous estimate of ODP. The bad estimate of ODP would then 
cause the state noise to be over or under-estimated depending on whether ODP was 
estimated high or low. In the case of an under-estimated ODP, this would lower the state 
noise Q k , which in turn would lower the Kalman gain. When the gain becomes small while 
the measurement still contains useful information, i.e., the crossing time of a satellite is 
changing due to an actual change in ODP, the filter is no longer functioning correctly and 
is said to diverge. Ultimately the filter should use the information contained in the 
measurement to correct the state, but when filter divergence occurs, the error covariance 
P k becomes small, indicating low uncertainty in the state, i.e., the current state is correct 
and the measurement erroneous. Thus, the new estimate of the state is not updated 
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properly to reflect the changes in ODP. In the case that ODP is over-estimated, having 
ODP in the state noise expression will drive Q k up and consequently both P k and the 
Kalman gam will increase. A filter with an artificially large gain is said to be reactive, with 
a high level of uncertainty in the state, i.e., the current state estimate is erroneous and the 
measurement correct, and therefore the filter reacts and corrects the state without 
accounting for measurement noise. 49 

Using a weighted averaging technique on ODP such as 


ODP, = (3QDP,_, + 20DP,_ 2 + QDP t J 

~6 


( 26 ) 


where ODP k represents the weighted average at time k, and ODP*.!, ODP*. 2 , etc., the 
estimated values for ODP at times k- 1, k- 2, ere., could reduce the lack of robustness or 
susceptibility to bad data. The use of a weighted ODP average in the state noise 
expression (25) would help prevent a single spurious estimate of ODP from making a 
large change to the state noise. On the contrary, if there were a real change in 
atmospheric density, ODP would tend to increase or decrease over consecutive 
measurements, which would shift the weighted average of ODP in a similar fashion. 

Measurement Noise Model 

To limit the number of variables that would need adjusting for filter “tuning” 
purposes, the measurement noise level was fixed, based on the typical variance offence 
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observations, and the ratio of the process noise to the measurement noise, Q*/R*, 
changed. The ratio Q*/R* is essentially the signal to noise ratio of the state variable with 
respect to the measurement, and is often used to tune a Kalman filter. 50 

Accountability of anomalous measurements is essential for successful Kalman 
filtering in a real environment, and ideally, bad data should be prevented from entering the 
filter in the first place. One method that could be used involves checking the 
“innovations” vector, [ z — Hx], for sudden jumps in amplitude or rate before a bad data 
point enters the filter. 51 

Process Flow 

The flowchart in Figure 24 depicts the Kalman filter process flow as applied to 
radar fence data. 
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Figure 24: Kalman filter process (low diagram 














The process begins in the upper left corner when an observation is obtained at some time 
T*. If mean orbital elements are available at some time T, between time T* and the 
previous observation T*. t , they are propagated in time to T*, otherwise the value of the 
mean elements at time T*.] are propagated to time T*. A fence-crossing time prediction is 
then made using the software PPT2. Next, the residual crossing error is computed 
(O - C)t . If the starting elements came from the DC process, the value of the residual is 
checked. If the crossing error is greater than +/- 2 seconds, the elements are discarded 
and the value of the mean elements at T*.i and x t _, are propagated to T k and a new 
prediction is calculated. If the crossing error is less than +/- 2 seconds, the DC elements 
and associated state, Xj , are kept. Subsequently, the state noise Q* is estimated using 

(24) or (25). The a priori information x~ and P“ is obtained by using O t _ u to 
propagate the state from x^jj to x“ using (16) and the error covariance from P*_, to Pa- 
using (17). The observation is then “processed” whereby the Kalman gain is calculated 
using (14) and both the state estimate x* and error covariance P* are updated using (13) 

and (15) respectfully. The value x * output from the filter is the value of n that should 
have been applied during the period from T*.j to T k to reduce the fence crossing error. 

The effects of n on the other elements, primarily t and n, must be accounted for by 

replacing the value of n or the state at T k .i with the new estimate x* and repropagating 
the elements to time T k . Updating the mean elements to T k concludes the recursive 
process until the satellite being tracked is observed again. 
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Objects Selected 


The Kalman filter process was evaluated using the three satellites listed in Table 5, 
which provides orbit and physical characteristics of the satellites. 


Table 5: Characteristics for addition al satellites 


Sat 

No. 

Int. 

No. 

Name 

Ha 

X 

Hp 

(km) 

i 

(deg.) 

RCS 

(m 2 ) 

Shape 

•3 

Size 

(m) 

Mass 

(kg) 

11848 

1980- 
051 A 

Meteor 1 
(30 th ) 

472 

X 

439 

74 

7.6 

Cyl. + 2 
vanes 

5.0 L ? x 
1.5 Dia. ? 

2200 

15326 

1984- 
104 A 

Cosmos 

1601 

1 458 

X 

436 

65.8 

14.8 

Cyl. 

4.0 L ? x 

2.0 Dia.? 

? 

16928 

1986- 
067 A 

Cosmos 

1776 

589 

X 

534 

97.7 

28.1 

Oct. 

Ellipsoid 

1.8 L ? x 
1.5 Dia. ? 

7 


The satellites were chosen based on the availability of radar fence observations and the 
level of residual crossing error. 


Satellite Residual Errors Prior to Filtering 

To determine if Kalman filtering can reduce the number of UCT’s, fence-crossing 
predictions were made for the three satellites during the first half of 1989 
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Day of Year, 1 989 

Figure 25: Residual crossing errors and ODP for Cosmos 1601 without filtering. 

without filtering. The top plot Figure 25 shows the residual crossing errors over time for 
Cosmos 1601 whereas the bottom plot is the corresponding values of ODP as calculated 
from DC of elements. As can be seen by examining the (O - C) plot, there are primarily 3 
periods when the crossing errors fall outside +/- 2 seconds for Cosmos 1601. The first of 
which corresponds to the large geomagnetic storm that occurred in March of 1989 when 
the geomagnetic index A p reached a value of 249, and where the crossing errors reached a 
maximum of nearly 12 seconds. Figure 26 shows a similar plot for Cosmos 1776. 
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Figure 26: Residual crossing errors and ODP for Cosmos 1776 without filtering. 


As before, the main disturbance occurred during March of 1989 along with several 
other periods where the satellite residual crossing errors fell outside the identification 
window limit of +/- 2 seconds. Notice that the ODP values were over-corrected during 
the geomagnetic storm on DOY 72, meaning that the ODP values derived from the DC 
solution were too large and caused fence crossing prediction times to be too early, and 
resulted in positive fence-crossing residuals outside the 2-second identification window. 
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Figure 27: Residual crossing errors and ODP for Meteor 1 without filtering. 

The top plot of Figure 27 shows the residual crossing errors for Meteor 1 . This 
satellite has only a few UCT’s suggesting that atmospheric drag was never large enough 
to seriously change the mean motion of Meteor 1 during the period of study. After 
examining the lower plot of ODP for Meteor 1 , it can be seen that the average magnitude 
of ODP for the period was approximately 3x10" revs/day 2 which is less than half of the 
nominal value recorded for the other two satellites. Referring back to Table 5, which lists 
RCS values and other characteristics for the three satellites, the RCS value for Meteor 1 is 
7.6 m‘, which is considerably smaller than the other two satellites. A lower RCS value 
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generally means lower surface area, which would tend to raise the ballistic coefficient of a 
satellite notwithstanding differences in mass. A higher P would lower the effect density 
and drag would have on orbital motion, which appears to be the case given the few 
number of UCT’s for Meteor 1 . 


Filtering Results 

Having established the residual crossing errors for the three satellites without 
recursive filtering, a direct comparison can now be made of the impact of filtering ODP 
during normal fence operations. A parametric study was performed using a wide range of 
state noise levels for each of the two different formulations of as delineated in equation 
(24) and (25). Figure 28 shows the effect of filtering the ODP of Cosmos 1601 using the 
state noise expression described in (25), where 20 percent of the previous estimate, 
ODP*.i, was used as an indication of the current state noise level of the atmosphere. 
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Figure 28: Effect of a one-state Kalman filter on Cosmos 1601 fence-crossing errors. 

The top plot demonstrates that the filter dramatically reduces the number of UCT’s with 
only a few points falling outside the +/- 2-second identification window. The bottom plot 
displays the values of ODP estimated by the Kalman filter for Cosmos 1601. Notice that 
after the geomagnetic storm on DOY 72, the filter estimate of ODP went negative, prior 
to recovering to a nominal value. Negative values of ODP are uncommon and could be 
due to a satellite undergoing maneuvers. In this case, the negative ODP is due to the filter 
over-correcting the state, meaning that too large a value for ODP was estimated during 
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the storm. When the geomagnetic disturbance subsided, ODP had to be reduced greatly 
to keep the fence-crossing errors to a minimum. Such wide changes to the state are 
necessary at both the onset and conclusion of solar storms. During these periods the filter 
state noise must increase for the filter to have a fast transient response to the storm 
enabling sharp corrections to ODP, but at the cost of adding more noise to the ODP 
estimate. 52 Because the Kalman filter is single state, fence-crossing etrors are essentially a 
function of ODP only, and therefore all the uncertainty is placed in ODP, rather than being 
distributed across additional state variables such as n and l . In actuality, some of the 
fence-crossing residuals are due to errors in n and l , and therefore a single state filter is 
susceptible to over-correction during severe solar storms. 

A number of trials were run using different levels of state noise for both the adaptive and 
non-adaptive state noise expressions for all three satellites. Table 6 and 
Table 7 summarizes the results. 


Table 6: Filter results for an adaptive state noise. 


Satellite 

Lo Filter 

in~^ — 1 

I — 


0.05 

0.10 

| 0.20 

0.40 

Cosmos 1601 

UCT's 

19 

16 

10 

4 

U ll 

9 


RMS 

1.26 

0.83 

0.67 

0.55 

0 65 

Cosmos 1776 

UCT’s 

18 

11 

9 

10 

r 12 


RMS 

1.27 

0.98 

0.89 

0.87 

1 17 

Meteor 1 

UCT's 

4 

3 

3 

6 

6 


RMS 

0.53 

0.52 

0.52 

0.55 

0.62 
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Table 6 shows the results of Kalman filtering on the number of UCT’s and on the root 
mean square (RMS) of the residual crossing error, using the adaptive state noise equation 
(25), with different values of a. The largest reduction in both UCT’s and RMS of the 
crossing error is shown by Cosmos 1601, followed by Cosmos 1776. Meteor 1 shows 
only a slight improvement over the no filter case. 


Table 7: Filter results for a non-adaptive state noise. 


1 Satellite 

No Filter 


) | 

l . 


1.00E-05 

1.00E-04 

2.00E-04 

4.00E-04 1 

Cosmos 1601 

UCT's 

19 

16 

6 

6 

l" i in 

12 


RMS 

1.26 

1.16 

0.59 

0.59 

0.76 

Cosmos 1776 

UCT’s 

18 

12 

9 

10 

13 


RMS 

1.27 

1.08 

0.90 

1.20 

1.27 

Meteor 1 

UCT’s 

4 

3 

6 

5 

4 


RMS 

0.53 

0.52 

0.94 

1.10 

0.53 


Table 7 shows the results of Kalman filtering using the non-adaptive state noise 
equation (24) with different values of y . In this case, the filter is able to achieve the same 
reduction in UCT’s and RMS of crossing error as in the adaptive state noise case, with the 
exception of Cosmos 1601, where the adaptive state noise expression lowered the number 
of UCT’s and RMS of crossing error further. 

These results suggest that the Kalman filter performs just as well using either the 
adaptive or non-adaptive state noise expressions for satellites with medium to low residual 
crossing errors. However, for satellites with large crossing errors, filtering using an 
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adaptive state noise appears to have greater success at lowering UCT’s and tightening the 
RMS of crossing errors. By having ODP in the state noise expression, Q* can increase 
during periods of severe solar storms, enabling the filter to make the larger, sharper 
corrections needed to ODP during such noisy periods, whereas during quiet atmospheric 
intervals, ODP will be smaller, making Q* contract, which will tend to make the state 
corrections smoother and less likely over-compensated due to noisy measurements. 
However, the non-adaptive Q* does not have the capability to expand or contract based on 
atmospheric conditions, and therefore the filter is more likely to have difficulty making 
transitions from noisy to quiet atmospheric conditions for satellites that are sensitive to 
large changes in density and drag such as Cosmos 1601 . 
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7. Conclusions and Recommendations 


These findings support the potential for estimating the orbital decay rate of 
satellites using an orbital decay parameter (ODP) model developed from fitting predicted 
density to actual ODP values in a least squares (LS) sense. For satellites in high eccentric 
orbits i.e., e > 0.2, or for satellites undergoing high decay rates, improvements can be 

I 

made in the ODP model by basing orbital decay estimates on the LS fit of pH 2 to satellite 
ODP rather than just density. Predicted ODP could then be used to calculate fence-crossing 
times and potentially reduce the number of uncorrelated targets (UCT’s) occurring at the radar 
fence. Such a system would be dependent on reasonable solar and geomagnetic predictions in 
order for atmospheric density to be properly estimated. Given these constraints though, 
estimating orbital decay rates using an ODP model would eliminate the phase lag that currently 
exists between the ODP calculated from EXT and changes in atmospheric density. This would 
prove to be critical during severe solar storms when up-to-date ODP values are necessary to 
successfully predict crossing-times and maintain the satellite catalog. 

Atmospheric Modeling 

Of the four atmospheric models compared, the Jacchia ’71 (J71) model performed 
the best as measured by the correlation between the density predicted by the atmospheric 
model and the actual ODP calculated from fence observations. Specifically, out of a dozen 
satellites examined, J7 1 showed the highest correlation for seven. The major factors behind 
the results were due to differences in how well each model estimated density fluctuations 
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due to variations in daily solar flux and geomagnetic activity. There was one particular 
geomagnetic storm that was greatly underestimated by all four models, clearly indicating 
the need for better modeling of the upper atmosphere. 

The ballistic coefficients of the satellites examined were determined using orbital 
decay theoty in conjunction with density estimated from an atmospheric model. Results 
indicated good agreement with theoretical values. Having the ability to determine 0 can 

help in identifying an unknown satellite or in detecting changes in mass or frontal area due 
to collisions or maneuvers. 

Kalman Filtering 

The use of a single state Kalman filter offers a potential improvement in 
maintaining the catalog of artificial satellites during large solar storms by reducing the 
number of UCT’s occurring and tightening residual RMS crossing ereors. LEO satellites 
with large surface-to-mass ratios whose ground-tracks extend into the geomagnetic pole 
regions of the Earth will be more affected by perturbations to atmospheric density as a 
result of solar activity and more likely to benefit from applying a Kalman filter to process 
their observations than satellites that are not as highly affected by drag such as Meteor I. 

An adaptive state noise expression capable of expanding the state noise during 
periods of severe atmospheric disturbance while contracting the state noise during nominal 
conditions, proved to reduce the number of UCT’s more effectively than a non-adaptive 
state noise for satellites susceptible to large crossing errors such as Cosmos 1601. The 
presence of ODP in the state noise formulation was the essential indicator of noise-level. 
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During periods of increasing solar or geomagnetic activity, ODP would grow as well as 
state noise, whereas immediately after a stormy period, ODP would decrease, bringing the 
state noise down as well. 

To lower the possibility of a filter diverging or becoming reactive as a result of a 
bad estimate entering the adaptive state noise expression, a weighted running-mean of 
ODP should be used. This would help prevent an over estimated value of ODP from 
opening up, i.e., increasing the state noise during otherwise quiet atmospheric conditions, 
which could lead to an over reactive filter, or in the case of a severely underestimated 
value of ODP, filter divergence. The non-adaptive state noise expression on the other 
hand, is not dependent on ODP and as a result, the filter is less vulnerable to divergence 
and reactivity problems. 

Regardless of which state noise expression is employed, the Kalman filter is still 
susceptible to bad data entering the process, which leads to poor estimates of ODP. This 
can be circumvented by checking each measurement and down-weighting, or increasing 
the measurement noise, for those observations that fall outside an expected range. This 
leads to a slower response from the filter to real changes in satellite ODP, but prevents the 
filter from making abrupt changes to ODP based on poor data. 

Implementing a Kalman filter system to help maintain the satellite catalog will most 
likely require individual filters for each satellite. Each filter would have to be tuned in 
order for the filter to respond correctly during real atmospheric events while ignoring 
noisy measurements. The process of tuning requires adjusting the state noise so that the 
filter can make corrections to the state that will minimize residual crossing errors without 
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causing the filter to diverge or become reactive. The degree to which a satellite is affected 
by drag would determine whether to use the adaptive or non-adaptive state noise 
expression. In some cases it may be possible to use the same state noise expression for 
satellites in similar orbits and with similar physical properties. 


Future Work 

The two research approaches investigated thus far suggest future work in the 
following areas: 

Real-Time Systems Operation 

Differential correction of elements needs to be performed in a real-time 
environment. This would allow testing the approaches of using atmospheric models 
and/or Kalman filters to improve estimates of ODP during periods of high solar activity, in 
preparation for incorporating these techniques into real-time operations. 

Multi-state Kalman Filtering 

The robustness of the Kalman filter technique needs to be improved. Using a 
multi-state Kalman filter, i.e., a filter having a state vector of multiple components such as 
£ ’ ", and ODP, to process fence-crossings as opposed to a single state, may provide a 
more accurate estimate of ODP. A one-state filter has only one degree of freedom, requiring 
the ODP output of the filter to solely eliminate residual crossing errors. This can lead to over 
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correction since there are other variables that could be responsible for the inaccuracy in fence 
crossing time, such as errors in £ and n. 

Hybrid Approach 

A hybrid approach to reducing UCT’s involving both the orbital decay parameter 
model and the Kalman filter process may prove beneficial. This technique would take 
advantage of the predictive capability of atmospheric models to estimate the initial state or 
a priori of a Kalman filter, rather than just guessing the initial value, which is common 
practice when initializing a Kalman filter. 
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Appendix A: Atmospheric Modeling 


To model the upper atmosphere, a number of assumptions are made. The 
atmosphere is generally assumed to be comprised of the following constituents up to 90 
km: molecular Nitrogen (N 2 ), molecular Oxygen (0 2 ), Argon (Ar), and Helium (He), 
homogeneously mixed with a fixed volume composition leading to a constant mean 
molecular weight M. f Above this altitude, dissociation of molecular oxygen due to 
extreme ultraviolet (EUV) absorption lowers the mean molecular weight. 53 The 
atmosphere is assumed to be in hydrostatic equilibrium yielding the relationship: 

dP = — gp dz ^ 27 ) 

where dP is the differential of pressure, p is the density, g is the height-dependent 

acceleration of gravity, and dz is the differential of geodetic height. In addition, the air is 
taken as an ideal gas with an equation of state 

T 


P = pR 


M 


(28) 


where R is the universal gas constant and T is the temperature. Upon substitution of (28) 
into (27), the barometric equation results 54 

Mg 


d/np = din 


T 


RT 


■ dz 


(29) 


After integrating and applying boundary conditions at 90 km, the density at altitude z, up 
to the homopause near 100 km, is given by (30). 


m = 

K T(z) 


exp 


'-fj^dz 


(30) 


90 


88 



Above this altitude, the departure from mixing and onset of molecular diffusion require the 
atmosphere to be modeled by the diffusion equation 55 


dn, 

n ; 


Mil 

RT 


dz - ( 1 + ctj )- 


dT 


(31) 


where n* represents the number density, or the number of molecules of the ith species per 
unit volume, Mi the molecular weight of the ith species, and 0Cj the thermal diffusion 
coefficient of the ith species. The solution to (31) is given by 56 


n j(z) = n,(100)| 

Total density is then found using 


T(100) 

T(z) 


Ml+ar,) 


exp 


R J T 

100 1 


(32) 


p = — M; 


(33) 


where N A is Avogadro’s number. These equations form the general underlying basis for 
determining atmospheric density, given a temperature profile T(z) and altitude z. 

There are a variety of temperature profiles in use in atmospheric models. Those 
typical of L.G. Jacchia, specifically J71, 57 begins at a boundary condition of Zo = 90 km, 
where the temperature starts at a fixed value of T 0 = 183° K, and has a gradient of 


G 


0 



(34) 


An inflection point occurs at a fixed height of z* = 125 km, above which the profile 
becomes asymptotic to a temperature T„, referred to as the exospheric temperature. The 
temperature at z x is given by 


f Minor constituents are ignored; Hydrogen is introduced between 150 - 500 km, depending on the model. 
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T x = a + bT_ + ce kT_ ( 

where a = 371.6678, b = 0.0518806, c = -294.3505, and k = -0.00216222. 

For temperatures in the region of 90 to 125 km, the temperature profile is given by the 


fourth- order polynomial 




subject to the following constraints 


when z = z. 


G. = 


when z = z. 


T - T 
1.9 — — -± 
z, - z n 


= 0 


The coefficients c n can then be solved for in terms of T x , and are given by 

. _ .„( T »-T„) 


c, = 1.9 


( Z x ~ z o) 


c 2 = 0 


c 3 = -1.7 


c 4 = -0.8 


(t, - T„) 
(z. - Zo)’ 

(t, -T 0 ) 

( 2 , - 2 o)' 


For altitudes above 125 km, temperature profiles are given by 
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T = T x + e tan' 


i 


(z- Zl )[l + tf(z-z x )*]j, 


where £ = — (T.-Tj, t? = 4.5 x 10 -6 and £, = 2.5. 


(38) 


91 



Appendix B: Model of Orbital Motion 

The model of orbital motion used to maintain the satellite catalog applies a 
technique similar to the method of variation of parameters, but does not require numerical 
integration. The model, PPT2 34 , employs an algorithm developed using Hamiltonian 
mechanics to propagate orbital elements in time by adding factors correcting for 
perturbations due to the non-sphericity of the Earth. PPT2 is based upon satellite orbital 
theory without drag as developed by Dirk Brouwer 58 , coupled with improvements made by 
R. Lyddane 59 that remove singularities at e = 0 and i = 0. In addition, drag is modeled 
using time derivatives of mean anomaly. What follows is a brief description of the 
Brouwer/Lyddane corrections to the elements, proceeded by the PPT2 drag model and a 
discussion of the satellite position prediction routine. 

Brouwer/Lyddane Model 

The Brouwer/Lyddane model 60 of orbital motion accounts for the non-sphericality 
of the Earth by using a spherical harmonic representation of the geopotential 

U = — + ~X~X^ m (si n /?)[C n m cos(mA) + S n m sin(mA)] (39) 

where (3 is the satellite latitude, X is the longitude, C„. m and S n . m are coefficients which 
depend on the mass distribution, P n m are the associated Legendre polynomials, p is the 
gravitational constant, R® is the equatorial radius of the Earth, and r is the magnitude of 
the position vector of the satellite. Equation (39) is often approximated by ignoring the 
longitudinal terms resulting in only a zonal harmonic approximation to the geopotential: 
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u 


(40) 



symmetric about the north-south axis of the Earth, i.e., they are not dependent on 
longitude. In addition, even-numbered harmonics are symmetric about the equatorial 
plane, whereas odd-numbered harmonics are anti-symmetric. 61 Due to the presence of 
inverse powers of r in the geopotential, equation (40) may be truncated. The Brouwer 
model takes advantage of this and uses an expansion of (40) in the first four zonals 


U = y + J 2 -^y-(l - 3sin 2 p) + J 3 -^y_(3sin (3 - 5sin 3 p) 

— ^—^-(3 — 30sin 2 P + 35sin 4 p) (41) 

8 r 

— — J 5 (l 5sin P - 70sin 3 p + 63sin 5 P ) 

8 r 

where J, = 0.4841605xl0' 3 V? J, = -O.QSOSSxlO' 6 -Jl 

J 4 = -0.55199x10-* J9 J 5 = - 0.65875x1 0“ 7 Vn 

The variation in the elements is separated into secular and periodic corrections. 

The secular corrections are functions of the even zonals, i.e., J 2 , J 4 , etc., eccentricity, semi- 

major axis, and inclination. Periodic corrections are broken down further into long and 

short period corrections. Short period corrections are functions of J 2 and the elements, 

while long period corrections are functions of all the zonal terms and elements. The 

secular and periodic corrections are listed, followed by the Brouwer/Lyddane algorithm 

for propagating orbital elements. 
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Secular Corrections 


Define tj = jl-e " 2 9 = co si". 


and introduce the following dimensionless variables 


r 2 = ^ r 3 =-4 t j , r3 


2a 


#3-3^© Yi - 


2 

8a' 


Y 5 =-- 


w 


y' — Zl. 
73 V 6 


T 


The secular corrections are represented by 


8,1 = |r ' 2 r](-l + 39 2 ) 


_ 7s 
7?’ 


7s = .,,0 


2 

+ 32 r ' 2 + 1677 + 25t] 2 + ( 30 ~ 96t? ~ 90r l 2 )0 2 + (105 + 14477 + 25t7 2 )0 4 ] 


15 


+— 74^ e/ ' 2 ( 3 - 3 °0 2 + 350 4 ) 


= — >"2 (— 1 + 5 a 2 ) 


+ 32 t " 35 + 2477 + 25772 + ( 9 ° “ 19277 ” 126r ? 2 ) 02 + (385 + 36077 + 45/7 2 )0 4 ] 


^"h r ' A [ 21 “ 9772 + (- 27O + 1267 r)0 2 +(385- I 8977 


2 l£)4 


2 

<5,Q = -3y' : 0 +- 72 2 [(-5 + 1277 + 977 2 )0 + (-35-36n-5T7 2 )0 3 ] 


■fr;(5-3rr)(3-70=)0 


94 



Long Period Corrections 


<5,e = e'V 


-72 

fi-iw’- 409 \1 

5r; r 


! - 50 2 _ 

12 y\ [ 


1 - 30 2 - 


80 4 


1 - 50 


7] 2 sin i" f 


4y' 
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2 240 


1 — 50 2 
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(\ , 
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>cos<a 
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+ e' 3 0 2 sin i*1 


. 320 2 
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Short Period Corrections 


8 2 z = a "y 2 


(-1 + 3 e 2 ) 


r' 3 77 3 ) 


*//3 

+ 3(l - 0 2 )— cos(2 co' + 2/') 


8 2 c = 


TYj 


2e' 
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= - ^-0 [6 (/' - t + e'sin/')- W2l] 


97 



Propagation Algorithm 


• Begin with the initial values of the mean elements a*, e', t cq*^ anj Q* 

Note: Double primes denote mean elements, which also may be obtained through 
suitable averaging of the osculating or instantaneous Kepler elements. 


• Calculate mean motion: n 0 = 

• Propagate mean anomaly, argument of perigee, and ascending node: 

*■' ~ £ o + n 0 t(l +SJ) 



co r = co' + n 0 t 
^ = + n 0 t 8 s Q 

where 5 S x represents the secular contribution to the element x. 
• Compute 8 , e, e'8, £, <fy\ and (sin/') 5, Q 

where 8 ]X is the long period correction to element x. 


(42) 

(43) 

(44) 


• Calculate 5iz = + 5,(0 + 8,n 

• Solve the system of equations 


E” - 

• e'sinE' = t" 

(45) 

1 

1 1 + e' 1 


tan-/' 
2 J 

= Vl-e' tan 2 E ' 

(46) 

rr 

a'(l-e" 2 ) 


r = 

(l +e'cos/') 

(47) 

for true anomaly and position. 



Calculate S?z = + 83(0 + 

82 Q 



98 




where ^x represents the short period corrections to element x. 

• Calculate z = £ + co + Q = £* + co" + Q." + 5iZ + biz 

• Compute z"8 £ = e / '<5,f + e'8 2 £ and Si = 8 { i + 

• Calculate 


1 




sin— i" J) Q = 


(sinz')(<5, £2 + 5 2 Q) 


2 cos— i' 
2 


<5,z 


• a = a' + 82 a 


( 48 ) 


• Solve the system below for e, £ , i , and £2. 

ecosC = (e" + 8 e)cos£" — t"8 £sin£ A 
esinf = (e" + 8 e)sin^" + e"8 £ cosf* 


/ 1 ^ 
sin— z 


cosQ = 


l ( 1 

sin-z + cos — z 


8 i 


cos£2* — ^sin^-z* , j8£2sin£2* (49) 


1 1 

sin £2 = 

1 

f 1 1 

1 „ " 


f 1 ) 

sin— i 

sin— i + 

cos —i* 

— <5 i 

sin£2* + 

sin — i" < 

2 J 


L 2 ! 

2 J 

2 J 


2 J 


• Determine the satellite position vector r using equations (45) - (47) with the newly 
found osculating elements. 
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Atmospheric Drag 


The satellite orbit propagator, PPT2 34 , does not contain an atmospheric model to 
handle drag effects on satellite orbital motion. The effects of drag, as well as other non- 
central forces such as solar radiation pressure, satellite thrusting, etc., that are not 

contained in PPT2, are represented by the time rate of change in mean motion as shown in 
equation (50), 

1 = f « +mt + ^ nt !+ | nf3 +- (50) 

where m = n 0 t(l +Sj), 

t is mean anomaly, is mean anomaly at epoch, i.e., t = 0, no is the mean motion at 

epoch, n is the time rate of change of mean motion, n is the time rate of change of n, 

and t is time. Normally, the n term is ignored, except for satellites with high decay rates. 

In addition, mean motion is updated using 

n = n 0 + nt (5I) 

where no is the mean motion at epoch. 

For a LEO satellite, n, or the orbital decay parameter (ODP), is primarily the 
result of atmospheric drag acting in a direction opposite to the satellite velocity vector. 
Changes in the mean motion of a satellite will effect the remaining orbital elements 
describing the orbit. Primarily, the effects are seen in the semi-major axis and eccentricity 
of the orbit as described by equation (52) and (53) below. 62 
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2 a • 

a n 

3 n 


2e , • 

e = - — (l-e 2 )n 
3n 


( 52 ) 


(53) 


Drag lowers the semi-major axis and eccentricity of a satellite orbit, thus tending to 
circularize the orbit. PPT2 uses (52) and (53) when propagating the mean elements to a 
specified time, as well as when making fence-crossing predictions. 


Fence-Crossing Time Prediction 

The capability is included in PPT2 to predict satellite fence crossing time for the 
purposes of preparing a chronological schedule of upcoming observations. When a 
satellite observation occurs, the actual crossing time is compared to the predicted, and if 
the observed minus the calculated (O-C), or residual error, is within two seconds, the 
satellite is considered to be identified. A diagram depicting the geometry of a typical fence 
station observation is presented next, followed by the method used in PPT2 to perform 
fence-crossing predictions. 

Station Observation Geometry 

Figure 29 illustrates the geometry of a typical radar fence observation, 


101 



z 

4 



Figure 29: Fence-crossing geometry 

where point A represents the ground station, R is the station vector, r is the satellite 
position vector, z is the local vertical, p is the slant range, AB is the ground range, BC 

is the vertical range, and 0 is the angle of elevation. Here the plane of the paper coincides 
with the fence plane. 

PPT2 Prediction Algorithm 

The capability to predict the geometry described in Figure 29 as well as fence 
crossing time is provided in PPT2. The condition that must be met signifying a satellite is 
crossing, i.e., in the radar fence plane is illustrated in Figure 30, 


102 



Polar Axis 



Figure 30: Satellite crossing radar fence-plane 
where a satellite with position vector r is crossing the radar fence depicted by the line 
perpendicular to the local horizontal and with unit normal vector n , and D is the fence- 
plane displacement from the center of the Earth. The crossing constraint can be expressed 
mathematically as 

n»r + D < S (54) 

meaning that a satellite is considered to be crossing the radar fence when the fence normal 
component of the satellite position vector is equal to the fence displacement from the 
center of the Earth within an error tolerance 5. 
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The method used to make fence-crossing predictions begins as follows. 

• Starting with the time of the elements or epoch, t = 0, update eccentricity using 

e' = e' + et 

• Compute the change in mean anomaly 


= E - e'sinE - £' 


where E 


'v/l-e* 2 sin / 
1 + e*cosf 


(56) 


Solve the cubic equation t 
starting with At = A £!m 


J = 


+ n 0 t + t 2 + -g- 1 3 for time using iteration. 


A£ - mAt - — At 2 + — At 3 
2 6 

At = At + — 

m 


Iterate while j J J > 10“ 4 up to 20 times, 
then t = t + At 


Update a", of, and Q" using 


(57) 


a = a 


a At 


dco” 

co = co" + ——At* 
d£ 


_ dor 
a = + 


• Apply periodic corrections to a, e, /, co, fi, and £ . 

• Solve Kepler’s equation (45), and update position using (47). 


Calculate 


dcy _ df m(l + ecos/) 2 

“■ = “< ~wBr 


(58) 


(59) 


104 



Compute the correction to true anomaly: 


A/ = 


d / (n* r + D) 
dt n • v 


(60) 


If A/ is less than a tolerance value, the prediction is finished. Otherwise /is adjusted 
by A/, and the procedure starts over again beginning with (55). 63 
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